#### Overlay density plots of BIOCLIM and microCLIM on the same plot space
#### Subset plots by bins of elevation
#### NB Ylims for plots created within plotting loop, then commented out, and read in a .csv file

### Load libraries

library('SDMTools')

### Establish directories

bc.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/BC6/'
mc.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/microCLIM/'
out.dir = '/home1/99/jc152199/climateimages/'
test.dir = '/home1/99/jc152199/testimages/rawdata/'
test2.dir = '/home1/99/jc152199/testimages/byelevation/'
test3.dir = '/home1/99/jc152199/testimages/byfpc/'
densitytest.dir= '/home1/99/jc152199/testimages/densitytest/'
setwd(densitytest.dir)

### List files in those directories

toimage = c(list.files(bc.dir,full.names=T)[1:4],list.files(mc.dir,full.names=T)[c(1,2,4,5)])

### Get a vector of names for the images

imagenames = gsub('.asc.gz','',c(list.files(bc.dir)[1:4],list.files(mc.dir)[c(1,2,4,5)]))

### Use a for loop to plot them all out to the same directory

for (i in c(1:8))

	{
	
	t.asc = read.asc.gz(paste(toimage[i]))
	
	png(paste(out.dir,imagenames[i],'.png',sep=''))
	
	image(t.asc)
	
	dev.off()
	
	cat(imagenames[i],'\n')
	
	}
	
# Close loop

### Read in some of the ASCII files on which the microclimate model was built

dem = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC/dem_WTplusbuffer_LatLong_WGS1984_250mres.asc')
dem.mod = dem
fpc = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC/fpc_WTplusbuffer_LatLong_WGS1984_250mres.asc')
fpc.mod = fpc

### Need to summarize dem and fpc into quartiles
### Establish bins for dem & fpc

dembreaks = c(0,400,800,1200,1600)
fpcbreaks = c(100,200)

### Reassign the value for each cell based on it's bin

dem.mod[which(dem.mod>=dembreaks[1] & dem.mod<dembreaks[2])]=1
dem.mod[which(dem.mod>=dembreaks[2] & dem.mod<dembreaks[3])]=2
dem.mod[which(dem.mod>=dembreaks[3] & dem.mod<dembreaks[4])]=3
dem.mod[which(dem.mod>=dembreaks[4] & dem.mod<dembreaks[5])]=4

fpc.mod[which(fpc.mod>=fpcbreaks[1] & fpc.mod<fpcbreaks[2])]=1
fpc.mod[which(fpc.mod>=fpcbreaks[2] & fpc.mod<fpcbreaks[3])]=2

### Read in the ylims object created in the loop

y.lims = read.csv('ylims.csv',header=T)

### First create histograms of each layer, with all 3 representations, using axes with the same scale

setwd(densitytest.dir)
#y.lims = NULL
#y.lims.out = NULL

for (i in c(1,4,5,6)) # Examine by bins of DEM

	{
	
	t.toimage = toimage[which(substr(toimage,nchar(toimage)-7,nchar(toimage)-7)==i)] ### These are the three ASCII's to compare
	
	t.imagenames = imagenames[which(substr(imagenames,nchar(imagenames),nchar(imagenames))==i)] ### Names for the three ASCII's
	
	bc.asc = read.asc.gz(paste(t.toimage[1],sep='')) ### Read in BC ASCII
	
	mc.asc = read.asc.gz(paste(t.toimage[2],sep='')) ### Read in microCLIM ASCII
	
	x.lims = range(range(bc.asc,na.rm=T),range(mc.asc,na.rm=T))
	x.lims = c(x.lims[1]-1,x.lims[2]+1)
	x.lims = round(x.lims,0)
	
	png(paste('Comparison Surface 0',i,'.png',sep=''),units='in',width=10,height=2.5,res=500)
	
	plot(bc.asc,type='n',main=paste('Climate Comparison - Surface 0',i,sep='')) # Set the plot space
	
	par(mfrow=c(1,4)) # Set the plot space to be a panel of 1 * 3
	
	h=1 # Position tracker
	
	for (bin in c(1:4))
		
		{
			
		bc.density = density(bc.asc[which(is.na(bc.asc)==F & dem.mod==bin)]) ### Create a density object from the BC ASCII
		
		mc.density = density(mc.asc[which(is.na(mc.asc)==F & dem.mod==bin)]) ### Create a density object from the microCLIM ASCII
		
		#y.lims = c(paste(i,sep=''),paste(bin,sep=''),y.lims,max(bc.density$y),max(mc.density$y)) # Run loop through once first without plotting, to create the ylims object
		
		#y.lims.out = rbind(y.lims.out,y.lims)
		
		#y.lims = NULL
		
		plot(bc.density$x,bc.density$y,xlim=c(x.lims[1],x.lims[2]),ylim = c(0,max(y.lims[which(y.lims$surface==i),c(3:4)])),ylab = 'Density',xlab = 'Temperature',main=paste('Elevation ',dembreaks[bin],' - ',dembreaks[bin+1],' m.a.s.l',sep=''),type='l',col='red') # Set the plot space and plot the first density object
		
		points(mc.density$x,mc.density$y,col='blue',type='l') # Overlay points of mc.density onto previous plot space
		
		cat('Elevation Bin ',bin,sep='','\n')
		
		}
			
	cat('Surface 0',i,sep='','\n')
			
	h=h+1
		
	dev.off()
	
	}
	
# Close loop
